C=====================================================================
C  Define pseudostates in Laguerre basis
C      f(n,l,r) = dnorm(n,l) * (2*alfl*r)**(l+1) * exp(-alfl*r) *
C                 * Laguerre(2*l+2;n-1;2*alfl*r)
C  by using formulas just for Laguerre polinomials.
C=====================================================================
      subroutine laguer (Z, L, alf, nmax, nload, enrg, c, h,
     >   dnorm, dnorm2, work, fac, nfac)
      implicit real*8 (a-h, o-z)
      implicit integer (i-n)
      dimension c(nload, nmax), h(nload, nmax), fac(0:nfac)
      dimension enrg(nmax), dnorm(nmax), dnorm2(nmax), work(nmax)
C
C INPUT:
C -----      
C  Z - charge of hydrogen like ion. 
C  L - orbital momentum.
C  alf - parameter of basis.   alfl = alf / (l+1).
C  nmax - number of pseudostates.
C  nload - load dimension of arries
C
C OUTPUT:
C ------
C  enrg - eigenvalies
C  H    - Hamiltonian's matrix.
C  C    - matrix of eigenvectors.
C  Work, dnorm2 - work space for the programe.
C
      L2 = 2 * L
      alfl = alf / dble(L+1)
      a2 = 2.0d0 * alfl
C
C     Define factorials as n!=dexp(fac(n))
C
      fac(0) = 0.0d0
      fact = 0.0d0
      tmp = 1.0d0
      ntmp = 2 * nmax + L2
      if (nfac .lt. ntmp) then
         print*, 'nfac is not enough to store array of factorials,'
         print*, 'nfac has to be more then   2 * NMAX + 2 * L =', ntmp,
     >      ' but  nfac =', nfac
         stop 'nfac has to be more then   2 * NMAX + 2 * L'
      end if 
      do n = 1, ntmp
         fact   = fact + dlog(dble(n))
         fac(n) = fact
      end do
C
C     Define normalization coeff. Dnorm 
C
      c1 = dsqrt(a2)
      c2 = 1.0d0 / a2
      do i = 1, nmax
         dnorm(i)  = c1 * dexp(0.5d0 * (fac(i - 1)  -  fac(L2 + 1 + i)))
         dnorm2(i) = c2 * dexp(fac(L2 + i)  -  fac(i - 1))
      end do 
C
C     Define matrix of Hamiltonian
C
      c2 = -a2 * a2 * 0.5d0
      do i = 1, nmax
         do j = 1, i
            sm2 = 0.0d0
            do jj = 1, j - 1
               do jjj = 1, min(i, jj)
                  sm2 = sm2 + dnorm2(jjj)
               end do
            end do
            sm1 = 0.0d0
            do jj = 1, min(i, j)
               sm1 = sm1 + dnorm2(jj)
            end do 
            diag = 0.d0
            if (i .eq. j)  diag = 0.25d0
            res = (dnorm(i) * dnorm(j) * ( (Z/alfl - dble(l+j)) * sm1 +
     >         sm2) + diag) * c2
            h(i, j) = res
            h(j, i) = res
         end do
      end do 
c$$$C         
c$$$C     The following part is used to calculate  < n I QHQ I m >
c$$$C     where   Q = I  -  I 1s >< 1s I
c$$$C
c$$$      nmax = nmax - 1
c$$$      print*, 'Warnig, LAGUER works for   QHQ - problem.', 
c$$$     >   'Q = I  -  I 1s >< 1s I'
c$$$      print*, 'alf, nmax, l:', alf, nmax, l
c$$$      do n = 1, nmax
c$$$         do m = 1, nmax
c$$$            h(n,m) = h(n + 1, m + 1)
c$$$         end do
c$$$      end do 
c$$$c
C  if matc = 0, then  only eigenvalues.
c
      matc = 1
      call rs(nload, nmax, h, enrg, matc, c, dnorm2, work, ierr)
      if (ierr .ne. 0)  then
         print*, 'Programe "RS" finished unnormaly, ierr =', ierr
         stop    'Programe "RS" finished unnormaly, ierr ='
      end if 
      return
      end
